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Abstract. A method to couple interparticle contact models with Stokesian dynamics (SD) is introduced 
to simulate colloidal aggregates under flow conditions. The contact model mimics both the elastic and 
plastic behavior of the cohesive connections between particles within clusters. Owing to this, clusters 
can maintain their structures under low stress while restructuring or even breakage may occur under 
sufficiently high stress conditions. SD is an efficient method to deal with the long-ranged and many-body 
nature of hydrodynamic interactions for low Reynolds number flows. By using such a coupled model, the 
restructuring of colloidal aggregates under shear flows with stepwise increasing shear rates was studied. 
Irreversible compaction occurs due to the increase of hydrodynamic stress on clusters. Results show that 
the greater part of the fractal clusters are compacted to rod-shaped packed structures, while the others 
show isotropic compaction. 

PACS. 82.70.Dd Colloids - 83.10.Rs Computer simulation of molecular and particle dynamics - 83.60.Rs 
Shear rate-dependent structure 



1 Introduction 

The mechanical properties of colloidal aggregates are of 
fundamental interest in science and technology. To classify 
particulate gels and to understand their rheological behav- 
iors is a key element. When attractive forces act among 
nano- or microscale particles, they form finite-sized clus- 
ters or a space-filling network. The latter shows a solid-like 
response to external stress, so it is regarded as a gel. In 
general, particulate gels are classified into two types ac- 
cording to the attraction strength between particles [ .]. 
If the attraction strength is sufficiently large, the par- 
ticle surfaces are deformed at the bonding point, caus- 
ing non-central forces [2]. In this case, Brownian forces 
neither cause debonding nor tangential displacements be- 
tween contacting particles. So, branched tenuous struc- 
tures formed in the aggregation process are maintained [■ !] . 
On the other hand, if the attraction strength is weaker, 
denser and multilinking local structures, such as tetrahe- 
dral connections, are seen [4]. This can be explained by 
tangential displacements due to Brownian forces. 

The tangential displacements between contacting par- 
ticles, i.e. sliding, rolling and torsion, play important roles 
in the structure formation and mechanical property of col- 
loidal aggregates. However, for nano or microscale parti- 
cles, it is not simple to characterize these interparticle 
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interactions. For example, characterization of rolling re- 
sistances requires elaborate experiments such as AFM [5] 
and optical tweezers [ ] . Though these direct observations 
have clearly proven the existence of tangential forces, the 
particles available for such measurements are restricted to 
certain sizes. This is why there is still no general method 
to fully characterize the contact forces in colloidal systems. 
An alternative approach of investigating colloidal aggre- 
gates is to develop simulation methods. In particular, phe- 
nomena at the mesoscopic level are expected to hold all 
necessary particle-scale information, and the comparison 
between simulations and experimental observations can be 
used for the characterization of contact forces. 

This work introduces a simulation method of cou- 
pling interparticle contact models and hydrodynamic in- 
teraction models. The contact model used in this work 
is similar to the one developed in granular physics [7- 
1 1], which is able to capture aggregates maintaining their 
structures under low stress while being restructured un- 
der high stress. The bond strength is assumed to be suf- 
ficiently larger than the thermal energy kBT, therefore 
Brownian forces are not considered. Instead, hydrody- 
namic stress induces restructuring of clusters. The hydro- 
dynamic interaction model employed here is Stokesian dy- 
namics (SD) [j j-j 4], which provides the relations between 
velocities of particles and the forces acting on them in the 
Stokes regime. The evaluation of the hydrodynamic inter- 
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actions is the most difficult and time consuming part of 
the simulation due to its long-ranged and many-body na- 
ture. SD is based on Faxen's law and multipole expansions 
to obtain the far-field mobility matrix, which can simulate 
particle disturbed flows with reasonable computational ef- 
fort. 

We apply this simulation method to investigate the re- 
structuring behavior of finite-sized tenuous clusters under 
flow conditions. Investigation of colloidal aggregates un- 
der flow conditions is a traditional problem in colloidal 
science. The original study on cluster sizes under shear 
flows dates back to almost a century ago [in], which con- 
sidered the cluster growth due to shear-induced collisions. 
To estimate equilibrium cluster sizes, one needs to know 
about the breakup mechanisms due to the hydrodynamic 
stress as well. Theoretical studies of this problem appeared 
after decades [16-18], and simulation studies of aggregate 
breakups have been appearing over recent years [l!)-2")]. 
Restructuring of clusters is an additional and challenging 
issue in this context, since it depends on details of the con- 
tact forces. In order to focus on restructuring behavior, a 
special situation is considered: the shear flow is increased 
in a stepwise, thus less abrupt, manner than in previous 
works. In this case, the clusters are hardly broken; in- 
stead they are reinforced by new bonds generated during 
the restructuring process. The time evolution of clusters 
is expected to reflect the nature of contact forces. Some 
characteristic restructuring behavior was observed in the 
following simulations. 

The contents of the paper are as follows: the used 
methods, the contact model and SD, are briefly described 
in sect. 2.1 and sect. 2.2. The coupling for the overdamped 
motion is formulated in sect. 2.3. The optimization for the 
dilute limit of aggregate suspensions is given in sect. 2.4. 
Approaches to study the problem are explained in sect. 2.5 
and sect. 2.6. After describing the parameters used for the 
simulations in sect. 3.1, the results are shown by consider- 
ing two main issues: (i) how does the imposed shear flow 
result in the compaction of aggregates? (sect. 3.2) (ii) what 
is the tendencies of the shape formation and orientation? 
(sect. 3.3) A discussion about the compaction in terms of 
consolidation is presented in 4.1, the observed tendencies 
in 4.2, and the hydrodynamic effect in 4.3. Finally, the 
outcome of the work is concluded in sect. 5. 

2 Method 

2.1 The contact model 

2.1.1 Model for the elasticity 

A simple contact model was employed to simulate spher- 
ical particles cohesively connected. The interaction be- 
tween two particles is described by a cohesive bond in- 
volving 4 types of degrees of freedom: normal (the center- 
to-center direction), sliding, bending^, and torsional dis- 

^ We use 'bending' instead of 'rolling', because they are 
equivalent except for a numerical factor, but 'bending' is more 
intuitive for dealing with deformations of colloidal aggregates. 



placements [2, 7-11, 26-32]. These relative displacements 
are expressed by using position vectors and vectors fixed 
at respective particles. So, a rotation of the frame of ref- 
erence does not affect the result, i.e. objectivity is sat- 
ified [11]. In this work, the Hookean force-displacement 
relationships are assumed for these degrees of freedom, 
which are characterized by the spring constants kjq, kg, 
fee, and k^. 

Normal displacement Let us suppose two spherical parti- 
cles i and j located at r^'^^ and r^^^. The center-to-center 
distance r^^'^") = jr'*' — r*^-'^] is changed by the normal el- 
ement of the acting force. A monodisperse system is con- 
sidered here, so the radius of particle is denoted by a. The 
force-displacement relation is given by 

fI^'-''' =A:N(r(''^)-2a)n(»^^"), (1) 

where n^'^'^'> = (r^-'-' — r('^)/r*^*'-'^ is the normal direction. 
Sliding displacement Sliding displacement is a tangential 
element of the relative displacement between particles 
with fixed orientations. In order to express the sliding dis- 
placement vector d^'^'^\ unit vectors fixed to each particle: 
and were introduced, called contact-point in- 

dicators in this paper (Fig. 1). Using the indicators, the 
positions of the original contact points are written as fol- 
lows: 

ri[^^) = a|(*'^'\ r^^i}^ = r^^'> +a^^^'''\ (2) 

When two particles get in contact, i.e. at the stress-free 
state, the contact points are the same ro^c"''' = ro.c.\ and 
the contact-point indicators are set to = n^"^'^^ [(a) 

in Fig. 1] . The sliding displacement vector is given by the 

projection of the deviation ro.c,^ — ro.c.^ onto the perpen- 
dicular bisector between the two particles: 

= a{zi^(''J') - (Z\^(*'^) • n(^'^))n(^'^')}, (3) 

where Zi^^''-?) = ^(J-*) — So, the force- displacement 

relation is given by 

fI^'^^ = ksS'^^l (4) 

Bending displacement Bending is a type of tangential dis- 
placement involving rotation, with the angle between the 
contact-point indicators quantifying this displacement. 
This angle is assumed to be small, so it can be approxi- 
mated by the norm of the vector product [^ (J-*' x (-|(»y ) ) | . 
Since it includes some torsional element, the bending angle 
vector cp^'^'^'> is obtained by subtracting the normal part: 

By using the bending angle vector, the moment-angle re- 
lation is given by 

M^^''^ ^ kBa^<p'-'''\ (6) 
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Torsional displacement Torsion is the rotational displace- 
ment around the normal vector n^^'^\ In order to express 
the torsional angle, another set of unit vectors fixed to 
each particle: rj^*'-'' and jy^^-'^, are introduced, called tor- 
sion indicators (Fig. 1). When two particles get in contact, 
i.e. at the stress- free state, they are set by choosing ones 
from the vectors being orthogonal to the normal vector: 
rji^;f) . n'^i-j) = q and t]'^^'^'^ ■ n'-^'^^ = 0, and parallel to 
each other rj'^*'-'^ — jy^^'*) [(a) in Fig. 1]. Since the torsional 
angle is also assumed to be small, it can be approximated 
by the norm of the vector product jr;'*'^' x rj'^'^^l. The tor- 
sional angle vector is defined as the normal element 
of their vector product: 

Qi^.j) = X jy^^'^*)) • n^'^3^}n^''^\ (7) 

By using the torsional angle vector, the moment-angle re- 
lation is given by 

M^'^') = fcTft'e^''^'^ (8) 




Fig. 1. The contact-point indicators and ^^^''\ and the 

torsion indicators r]^''^^ and tj'--'''^ are illustrated for the stress- 
free state (a) and stressed state (b), respectively. The normal 
vector n^''-'' always indicates the center-to-center direction. 



configuration to release the potential energy stored in the 
tangential springs. 

The supportable strength for a cohesive bond depends 
on the direction of the acting forces or moments. So, the 
breakableness is characterized by two critical forces and 
two critical moments: i^No -Fsc -^Bc, and Mtc- In general, 
all components of the bond are stressed simultaneously. 
This is why a criterion of breakage can be introduced by 
a destruction function C{Fn, Fs, Mb, Mt), whose positive 
value indicates breakage. Here, a simple energy-like func- 
tion is used [ ■]: 

C = ^(^N)|f + 11 + ^ + ^-1, (10) 

^Nc -^Sc ^^'■'Bc ^''■'Tc 

where t?(i^N) is Heaviside function i?(Fn) ~ 1 for Fn > 
and i?(Fn) = for Fn < 0. 

According to the intensive studies by Dominik and 
Tielens [ ], the critical normal and sliding forces, Fnc 
and Fgc , are much larger than the corresponding forces of 
the critical bending and torsional moments, M^c/o, and 
Mtc/o- For a typical case, the ratio can be the order of 
10^. Thus, this work focuses on the effects for the bending 
and torsional breakups and excludes the separation and 
sliding breakups. Besides, the direct measurements of the 
critical bending moment have been reported [5, (i], while 
no direct measurement is available for the critical torsional 
moment. For simplicity the same strength for the bending 
and torsional moments is assumed here. In short, a spe- 
cial case of the bond breakableness written by Fnc — oo 
and — >■ oo and Mbc = -Mtc = -^c is considered. The 
strength of bond is then given with one parameter Mc by 

C-^'^-i. (11) 



Thus, the forces and moments on the contact point 
between particle i and j are related to the corresponding 
displacements. The force and torque acting on the particle 
i from the contacting particles j are given by their sums: 

^ (9) 

3 

The suffix P indicates the particle contact interactions in 
contrast to the hydrodynamic interactions. 

2.1.2 Model for the plasticity 

In the contact model, the potential energy is stored in 
the introduced bonds as long as the stresses acting on 
the bonds are small. When stresses become larger than a 
certain threshold, the bond breaks and the stored energy 
is dissipated. If particles are still in contact, the contact- 
point indicators and torsion indicators are reset with the 



2.1.3 Model for the new connection 

As of now, interactions between contacting particles are 
defined, but no assumption about particles being initially 
farther apart has been made. We consider a short-range 
cohesive interaction in this work. As a simple case, it is 
assumed that no interaction acts between remote parti- 
cles (r'^*'-'^ > 2a). If, however, two particles approach each 
other and get into contact, i.e. r^''-'' = 2a, they start to 
interact with each other, which is modeled by the genera- 
tion of a cohesive bond. 

2.2 Hydrodynamic interaction (Stokesian dynamics) 

Stokesian dynamics (SD) is employed for evaluating the 
hydrodynamic interactions [12-14]. Here, simple shear 
flows u°°{r) = z-yCj, are considered to apply, where 7 
is the shear rate. The force-torque-stresslet (FTS) version 
of SD is required to solve the flow conditions. By using 
the translational velocity C/°°, vorticity fl°° , and rate-of- 
strain E°°, the flow field M°°(r) is expressed as follows: 

w°°(r) = (7°° + X r -f E°°r. (12) 
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with the following nonzero elements: H!^ = j/2 and 
= = j/2. The hydrodynamic interactions act- 
ing on a particle i, i.e. the drag force f'^ , torque , 

and stresslet , are given as linear combinations of the 
relative velocities from the imposed flow: the translational 
and rotational velocities iJ^J) - u°°{r'-j'^) and f}'^^'^ - f2°° , 
of all particles {j = 1, . . . , N) and the rate of strain —E°°. 
The linear combinations for all particles are expressed as 
a matrix form 



(13) 



where the vectors involve llA^ elements for all particles, 
and the matrix R is the so-called grand resistance matrix-^. 

It must be noted that the lubrication correction of SD 
is not applied in this work. For suspensions where the in- 
terparticle interaction is absent, the lubrication forces play 
essential role for near contact particles [33, 34]. On the 
other hand, for rigid clusters, i.e. if the relative velocities 
between particles are zero due to strong cohesive forces, 
the lubrication correction has no contribution. Thus, the 
lubrication correction to the mobility matrix can safely be 
omitted [35-37]. Though, the relative velocities between 
particles are not perfectly zero in this work, near rigid- 
motion of clusters is investigated by only gradually in- 
creasing the shear rate. Due to the resulting small relative 
velocities of the primary particles, the lubrication correc- 
tion is expected to be less important. The neglect of lubri- 
cation forces is also a necessity for the computational ap- 
proach presented later. The investigated simulation times 
are only accessible with reasonable computational effort, 
if the time-scales of the long-range hydrodynamic forces 
and short-range contact forces can be separated. This sep- 
arability allows the reuse of the mobility matrix for sev- 
eral time steps which significantly enhances computational 
performance. This would unfortunately not be the case if 
lubrication forces would be considered. 



2.3 Overdamped motion 

To simulate the time evolution of particles with contact 
models, configurations of spherical particles are described 
by not only their central positions r(*^(t) (i = 1, . . . ,iV), 
but also by their orientations. The orientation of a par- 
ticle i is expressed by using a quaternion ^^'•'(t). If we 
set q^'^^O) = 1 at the initial time, the quaternion q^'^\t) 
means the rotation from the initial orientation. The ro- 
tation of a vector ^ fixed to the particle is written as 
|(t) = g«(t)^(0){gW(t)}-i. 

When the contact forces are strong enough, Brownian 
forces are negligible. But, the hydrodynamic forces depend 
on the imposed shear flow. Therefore, only contact and hy- 
drodynamic forces acting on the particles were considered. 



^ Since both the stresslet and rate-of-strain tensors are sym- 
metric and traceless, the five independent elements are denoted 



In general, the particles follow the Newton's equations of 
motion: 



dU ^ ^ dn 

m— = Fp + Fh, I— 
dt dt 



(14) 



where m and / are the mass and moment of inertia of 
the particles, respectively. The velocities C7, angular ve- 
locities fl, forces F and torques T include N vectors for 
all particles. 

For colloidal systems, the inertia of particles are neg- 
ligibly small compared to the hydrodynamic forces. By 
neglecting the inertia terms, the equations of motion (14) 
are approximated by the force- and torque-balance equa- 
tions: 

Fp + Fh « 0, Tp + Th « 0. (15) 

Systems following these balance equations are called over- 
damped. In order to solve the overdamped motion with 
SD, the mobility form: 



(16) 



is used instead of the resistance form (13). The numeri- 
cal library developed by Ichiki [ ] was used to obtain the 
mobility matrix M. By combining (15) and (16), the ve- 
locities of the particles {U, f2) are given by the functions 
of contact interactions {Fp,Tp): 





U{t) = U{Fp,Tp), n{t) - r2(Fp, Tp) 



(17) 



Once their velocities are determined, the time evolution of 
the particles are given by integrating the time derivative 
relations: 



d< 



dg 



nii) 



At 



(18) 



where the matrix Q^'^ is constructed by the elements of 
the angular velocity ^7*^*-' as follows: 



/ 







('0 
y 

(0 



-a 



(0 



-a 





a 



a 



(19) 



/ 



as vector forms, such as S = {Sxx, Sxy, Sx 



, Syz, Syy). 



Since the overdamped motions with the simplified con- 
tact model and approximated hydrodynamics are consid- 
ered, the accuracy of the numerical integration has no pri- 
mary importance. Therefore, the explicit Euler method 
was used to integrate the differential equations (18) with 
a discretized time step St. 



2.4 Reusing the mobility matrix for deforming clusters 

The bottle neck to simulate the time evolution is the cal- 
culation of the mobility matrix M in (16) in each time 
step. Since the contact forces are changed by short dis- 
placements of particles, the time step St needs to be set 
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small enough, causing a large calculational effort. This is 
why a way to reduce the computational effort has to be 
introduced. 

The mobility matrix M depends only on the positions 
of particles. If the relative positions of particles within 
an isolated cluster remain unchanged, the hydrodynamic 
interactions under any flow written as in (12) can be eval- 
uated with a single mobility matrix. Though clusters are 
not rigid in this work, up to a certain degree the deformed 
structure can be considered the same for the hydrody- 
namic interactions. As long as the deformation is negligi- 
ble in this sense, a mobility matrix may be reused repeat- 
edly. 

In order to evaluate the motion of an isolated cluster, 
one can take the centcr-of-mass of a cluster as the origin 
of the coordinate without loss of generality. Let us sup- 
pose that the deformation of the structure of the cluster 
for a time interval At is negligible. In this case, the time 
evolution of the particles from t to t' = t + At can be 
approximated by 

r«(t')«Rt^t'r-«(t), (20) 

where Rt^t' is a rotation matrix. The hydrodynamic 
interaction at the time t' , i.e. the relations between 

{F'i^\t'),T'^\t')) and (U^'^f), f2^'^t')) can be obtained 
by using the mobility matrix at the time t as follows: 

/AU{t')\ ( Fn{t') \ 

Af2{t') = -M{t) T^it') (21) 

where one has the following relations: 





\t') 


= R-4,,F«(t'), 


^ H 


\t') 






'it') 


p — 1 noo o 



and 

C/(*)(t') = Rt-,t'AU^"^ + U°°{r^\t')), ^^^^ 
r2W(t') = Rt-,t' Af}^''^ + f2°^\ 

Now one needs to determine the rotation matrix Rt^t' 
for the cluster, which is deformed during the actual time 
evolution. For a trial rotation matrix R, the positions of 
particles r*^*^ (t) are transformed to 

s(') = RrW(i). (24) 

The optimal rotation matrix Ropt should minimize the 
differences between the actual positions r^^\t') and the 
transformed positions s^*-*. One can take the following ob- 
jective function to be minimized: 



The gradient descent method is employed to find the opti- 
mal rotation matrix Ropt- Thus, the rotation matrix Rt^t' 
can be determined: Rt-ft' — Ropt- 

The objective function with the optimal rotation 
D(Rt_j.t') represents the degree of the deformation. If the 
deformation of the cluster becomes larger than a thresh- 
old: Z)(Rt_j.(') > I?maxi the mobility matrix needs to be 
updated. 

2.5 Introduction of a dimensionless shear rate 

The behavior of a cluster formed by strong cohesion un- 
der a strong flow is equivalent to the case of weak cohe- 
sion and a weak flow. In order to reduce the redundancy, 
a dimensionless variable, the ratio between hydrodynamic 
interactions and contact forces, is introduced. The cohe- 
sive force gives the typical force of the simulation Fq] the 
critical force for bending and torsional breakages is taken 
for that: Fq = Mc/a, because they play an important role 
for the restructuring of tenuous clusters. Since hydrody- 
namic interactions are proportional to the shear rate 7 in 
the Stokes regime, the dimensionless shear rate can be de- 
fined hy r = 67r?7oa^7/Fo = 67r7yoa^7/Mc, which indicates 
the flow strength for the contact force. In this work, the 
shear-rate dependence is discussed in terms of this dimen- 
sionless variable F . 

2.6 Stepwise increase of shear rates 

In general, three types of behaviors are expected for a 
cluster in shear flows: 

Rigid body rotation When the hydrodynamic stress is suf- 
ficiently weak, the cluster rotates without changing its 
structure. 

Restructuring When the hydrodynamic stress slightly ex- 
ceeds the strength of the cluster, the cluster is restruc- 
tured. Newly generated cohesive bonds during the restruc- 
turing may reinforce the cluster. If the strength of the 
cluster exceeds the hydrodynamic stress, it turns to the 
'rigid body rotation' regime. 

Breakup When the hydrodynamic stress is much stronger 
than the strength of the cluster, the cluster is significantly 
elongated and may be broken up into smaller pieces. 

In other simulation studies [19-21, 23-25, 39, 40], the 
shear flow is abruptly applied as a step function of time. In 
that case, the restructuring plays a limited role in a certain 
range of shear rates. This change of shear rate in a single 
step is a simple but very special case in terms of shear his- 
tory. If the flow strength is less abruptly increased, restruc- 
turing may reinforce the cluster before reaching higher 
shear rates. 

In this work, the focus is placed on such restructur- 
ing and consolidation aspects. So, the flow is turned on 
less abruptly. In order to plot the intermediate states of 
clusters by shear rates, the shear rate is increased in a 
stepwise manner. The fc-th shear rate is given by 

Fk = ri{rm../r,y'~'^^^''—-'\ (26) 
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where Fi is the initial shear rate, /Inax the final shear rate 
and fcmax the number of steps, and it is kept for the time 
period t^. resulting in an equivalent total shear strain F* , 
I.e. tl = F*/rk (see Fig. 2). 



O 



r = / f (0 At 



Fig. 2. The shear rate F is increased in a stepwise manner. 
The horizontal axis shows the total shear strain. 



regime with the lower shear rates and the sufficient com- 
paction with the higher shear rates. 

For reusing the mobility matrix for deformed clusters 
(see sect. 2.4), the threshold Dmax = O.Ola^ was given, 
which is small enough to evaluate the drag forces in ac- 
ceptable precision for our purpose. The actual update 
numbers of the mobility matrix during one cluster rotation 
are given in Fig. 4. These numbers are much less than the 
time steps to integrate the equations of motion, but these 
updates sufficiently reflect the long-range hydrodynamic 
interactions acting on deforming clusters. 



Table 1. Parameters of the imposed flows. 





Symbol 


SD 


FDA 


Initial shear rate 


A 


0.003 


0.001 


Final shear rate 


r 

^ max 


15.9 


10 


Number of steps 


k 

"'max 


28 


30 


Total shear strain 




20 


20 


on time interval 



3 Results 

3.1 Parameters for the simulation 

Fractal clusters generated by the reaction limited hier- 
archical cluster-cluster aggregation (CCA) were used as 
an initial configuration [41, 42]. The fractal dimension is 
df « 2. It is worth noting that such generated clusters have 
no loop structure. In previous works [37, 4.3], the hydrody- 
namic behavior of various sizes of the same CCA clusters 
have been examined by assuming rigid structures. Here, 
the restructuring behavior of small clusters with iV = 64 
was investigated. 

For randomly structured clusters, one needs to eval- 
uate a sufficient number of samples to study any gener- 
alizable behavior, therefore 50 independent clusters were 
simulated under the same conditions. A random selection 
of the initial clusters is shown as projections on x-z and 
x-y planes in Fig. 3 (a) and (b). 

The required parameters for the contact model (see 
sect. 2.1) are only the ratios between the spring constants 
of different modes and the critical moment. For the spring 
constants, the same value was set for the bending and 
torsional modes, and 10 times larger for normal and sliding 
modes: 

fcT = fcB, fcN = fcs = lOfcB. (27) 

The critical moment Mc for bending and torsional springs 
in (11) was set to the value which gives 1% of the particle's 
radius for the critical displacements. 

The used parameters for the imposed shear flows (see 
sect. 2.6) are presented by Table 1. In order to distinguish 
the hydrodynamic effect with Stokesian dynamics (SD), 
the free-draining approximation (FDA) was also used as 
the reference. For both methods, the ranges of the shear- 
rate changes were chosen to see the rigid-body rotation 
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Fig. 4. The mobility matrix is updated when the cluster de- 
formation exceeds the threshold. The average numbers of the 
updates during one cluster rotation are shown. The error bars 
exhibit the standard deviations over 50 independent simula- 
tions. 



3.2 Compaction 

First, the relation between the compaction and the flow 
strength was considered. The radius of gyration 

R>^j:ir'''-ro)\ (28) 

1=1 

where rp is the center of mass of the cluster, approxi- 
mately represents the hydrodynamic radius of the fractal 
clusters [37, 44, 45]. Indeed, the radius of gyration has 
been commonly used to quantify the size of random struc- 
tured colloidal aggregates. Fig. 5 (a) shows the shear-rate 
dependence of the radius of gyration, where final values 
at each shear-rate step were sampled. The averages and 
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~, N = 64) are shown by x-z and y-z projections (a) and (b), 
3 shown by x-z and y-z projections (c) and (d). 



standard deviations were taken over 50 independent sim- 
ulations. However, this quantity is not optimal to address 
the compaction behavior because results for compacted 
clusters depend on their shapes. 

The volume fraction is an alternative to quantify the 
compaction as it takes into account cluster shapes. As seen 
in Fig. 3 (c) and (d) , some of compacted clusters exhibit 
elongated shapes. Though the definition of volume frac- 
tion is not simple for isolated clusters, a rough estimation 
was used here. An arbitrary shaped cluster can be trans- 
lated into an ellipsoid having the equivalent moments-of- 
inertia. We take the ratio between the total volume of 
particles and the volume of the equivalent ellipsoid as the 
effective volume fraction (j)(.s (see Appendix A). The shear- 
rate dependence of the effective volume fraction is shown 
in Fig. 5 (b). The standard deviations are reduced in com- 
parison with the plot for the radius of gyration. So, the 
effective volume fraction is used as a measure for the com- 
paction behavior. 

For lower shear rates, both i?g and (j)f,s are almost 
constant, which indicates the 'rigid body rotation' regime 
described in sect. 2.6. By increasing the shear rate, the 



compaction starts at a certain point. This critical shear 
rate is denoted by I^c- For reference, the results by using 
the free-draining approximation (FDA) are also plotted by 
the triangle marks (A) in Fig. 5. The critical shear rate 
by FDA is much smaller, because the imposed flows are 
not disturbed by the cluster in FDA. Though the cluster 
size is small {N = 64), the ratio of the critical shear rates 

for the two methods is significant r^c^^ / Trc^''^^ ~ 4.2. 

Under stepwise increasing shear flows, clusters are 
monotonously compacted. The shear-rate dependence of 
the effective volume fraction turned out to be rather small 
so that the maximum compaction (d(/)eff/d/^ — 0) was 
not observed within the simulated range of the shear 
rates. The higher shear rate required for further com- 
paction may violate the model assumptions, such as Stokes 
regime for the hydrodynamics (sect. 2.2) and the condi- 
tions for the overdamped motion (sect. 2.3). The max- 
imum compaction in the SD simulations is (j)cg ~ 0.49 
with 7^(^^) w 15.9. The equivalent compaction by using 
FDA requires a shear rate of ^(fda) ^ ^q.G. It is worth 
noting that the range of shear rates for this compaction 
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Fig. 5. The compaction behavior is seen by the shear-rate 
dependence of the radius of gyration _Rg (a), and the effective 
volume fraction ^cff (b). The final values at the each shear-rate 
step were sampled, and the averages and standard deviations 
were taken over 50 independent simulations. The results with 
SD and FDA are shown by circles (Q) and triangles (A), re- 
spectively. 

is quite wide, i.e. f^^^^cj) = 0.49) /fjc^'^ « 1.1 x 10^. A 
discussion about this point will be given later (sect. 4.1). 

3.3 Shape and orientation during compaction 

The simulation results showed some characteristic ways 
of the clusters' compaction. In Fig. 3, the initial config- 
urations, (a) and (b), and the corresponding compacted 
clusters, (c) and (d), are displayed. It can be noticed that 
the initial clusters have a variety of shapes and orienta- 
tions, while the compacted clusters can be roughly clas- 
sified to two types: rod-shaped clusters elongating to y 
axis and round-shaped clusters. This section deals with 
the transitions of the shapes and orientations during the 
compaction. 

In order to quantify the slenderness of clusters, the 
aspect ratio of the equivalent ellipsoid r is used (see Ap- 
pendix A). The orientations can be quantified by the tilt- 
ing angle 0, which is the angle between the principal axis 
rii of the principal moment of inertia with the smallest 
value and y axis: 



= arccos(e„ • ni) 



(29) 



Now, the transitions of the shape and orientation of the 
distinct simulations can be represented by the aspect ratio 
vs. tilting angle {r-0) scatter plot. The four stages of the 



compaction, i.e. the initial (a), two intermediate (b) and 
(c), and final (d) configurations, are shown in Fig. 6. 

As mentioned above, the clusters seem to be distin- 
guished by two types after the compaction. In order to 
follow the formation processes, they are separated into 
two groups by introducing a selection criterion. The clus- 
ters having aspect ratios r > 1.56 at the results with 
r = 15.9 are classified as rod-shaped clusters and the oth- 
ers as round-shaped clusters. The threshold value r — 1.56 
was chosen by considering the distribution of the orien- 
tations, i.e. clusters having smaller aspect ratio than this 
value no longer show an orientation tendency. For the sim- 
ulated clusters, 74% of them are classified as rod-shaped. 
In Fig. 6, the circle (Q) and square (□) plots indicate the 
clusters ending up as the rod-shaped and round-shaped 
clusters, respectively. 

The representation for the initial clusters [Fig. 6 (a)] 
shows the following: The CCA clusters originally have 
anisotropic structures, so that the aspect ratios are dis- 
tributed between 1.5 and 3.5. Their orientations, on the 
other hand, are uniformly distributed. In this represen- 
tation, the two groups appears to be randomly mixed. 
Thus, at the moment we cannot predict which mode of 
compaction from the initial configuration. 

The compaction processes can be followed by seeing 
the translations of the plots from (a) to (d) in Fig. 6. 
Besides, in order to see the shear-rate dependence, the 
averages and standard deviations of the respective quan- 
tities by each group are shown in Fig. 7. For the effective 
volume fractions (^eff, the difference between two groups 
is not significant [Fig. 7 (a)]. For the clusters ending up as 
rod-shaped clusters, the distribution of the aspect ratio r 
is slightly compressed to smaller values at the beginning, 
and it is shifted to larger values after F ^ 0.14 [Fig. 7 (b)]. 
The distribution of the orientation shows a more system- 
atic change, i.e. the reorientation seems to correlate with 
the compaction [Fig. 7 (c)]. For the clusters ending up as 
the round-shaped clusters, the aspect ratio is decreased 
for r < 0.2, and stays constant after that [Fig. 7 (b)]. 

For reference, the results with FDA are also shown in 
Fig. 8. The observed tendency seen in Fig. 6 (d) is less 
clear in the equivalent compaction with FDA [Fig. 8 (a)] . 
For the same selection criteria, 42% of the clusters are 
classified as the rod-shaped group. The averages of the 
aspect ratio are taken over each group, and the shear-rate 
dependence is presented in Fig. 8 (b). For the result of 
FDA, the aspect-ratio of the rod-shaped clusters does not 
increase as in the SD simulations (shown by the dashed 
line). Instead it remains almost constant. Thus, by com- 
paring simulations of SD and FDA, it turns out that the 
features seen with SD are less pronounced with FDA. 



4 Discussion 
4.1 Consolidation 

As seen in 3.2, the shear rates for the compaction of clus- 
ters range over multiple orders of magnitude, i.e. initial 
fractal clusters are fragile while compacted clusters are 
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Fig. 6. The scatter plots of the aspect ratio vs. tilting angle {r-0) representations for the initial configurations (a), for the 
results at two intermediate shear rates F = 0.99 x 10~^ (b) and F — 0.91 (c), and for the final configuration with F — 1.59 x 10 
(d), show the restructuring processes of 50 simulations. The circle (Q) and square (□) plots indicate the groups of rod-shaped 
and round-shaped, respectively. For (b)-(d), the final values at the each shear-rate step were sampled. 
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Fig. 7. The shear-rate dependence of the effective volume fraction 0cft (a), aspect ratio r (b), and tilting angle (c), are 
shown by the two groups: rod-shaped and round-shaped, respectively. The error bars indicate the standard deviations over the 
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Fig. 8. (a) The aspect ratio vs. tilting angle (r-O) represen- 
tations for the compacted clusters {F = 10) by using FDA 
are shown. The circle (Q) and square (□) plots indicate the 
clusters ending up as the rod-shaped and round-shaped clus- 
ters, respectively, (b) The aspect ratios r are averaged over the 
members of the respective groups, which show the shear-rate 
dependence. The results by using SD are also shown by the 
dashed lines to compare with FDA. 



more and more robust to imposed flows. In our simula- 
tion with Stokesian dynamics (SD), the highest shear rate 
(resulting in (^cff w 0.49) is about 10'^ times larger than 
the critical shear rate /Ic where (j)ctf ~ 0.12. This may be 
explained by the following two non-linear effects: 



(i) The smaller the hydrodynamic radius is, the weaker 
the hydrodynamic stress acting on the cluster becomes 

(ii) The higher the number of newly generated loops within 
the cluster, the more the cluster resists restructuring. 

In order to see effect (i), the shear-rate dependence 
of the stresslet acting on a cluster was evaluated. The 
stresslet acting on a cluster is composed as follows [-Mi, '.^7]: 



«W ® F« + (iW (g)FW)T zW.fW 



where i^*^ = r*^'^ — tq- The stresslet S^'* for individual par- 
ticles i has already been calculated in (16). This stresslet 
Sci indicates the contribution of a single cluster to the 
bulk stress of a sheared suspension. For dilute suspensions 
of rigid spheres, the stresslet is proportional to the shear 
rate: (Ssph)a;z = (20/3)7r?7oa'^7- This is why the effect of 
restructuring on the hydrodynamic stress appears in the 
ratio between {Sci)xz and F. With the shrinkage of clus- 
ters, the efficiency is decreased [Fig. 9 (a)]. However, this 
decrease is not large enough to explain the wideness of the 
shear-rate range. 

The volume-fraction dependence of the stresslet was 
also plotted to see effect (ii) [Fig. 9 (b)]. If the total 
strain on time interval F* is infinitely large, the com- 
paction of cluster at each shear rate is expected to be 
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settled after some restructuring, and the further com- 
paction requires a higher hydrodynamic stress. This is 
why this plot approximately represents the compactive 
strength (yield stress) as a function of the volume fraction. 
Though the hydrodynamic stress given by the stresslet 
is not the same quantity as the mechanical stress, one 
may notice some similarity with volume-fraction depen- 
dent compressive yield stress Py{(t>) of space-filling col- 
loidal aggregate networks [46, 47]. The power-law behav- 
ior seen in the intermediate range of Fig. 9 (b) shows the 
large exponent 4.5, which was evaluated by fitting the av- 
eraged data within 0.13 < (^eff < 0.34. According to the 
cited work [47], large exponents of power-law relations 
were explained as a result of network densification due 
to irreversible restructuring. If the hydrodynamic stresses 
induce deformation of clusters, the same consolidation ef- 
fect is expected by the rule of new bond generations given 
in sect. 2.1.3. Thus, it was confirmed that the effect (ii) is 
responsible to explain the wideness of the shear-rate range 
for the compaction. 




Fig. 9. (a) The shear-rate dependence of the ratio between 
the stresslet (Sci)i:z and the shear rate F are shown, (b) The 
volume-fraction dependence of the stresslets acting on clusters 
are plotted. For both the plots, the final values at the each 
shear-rate step were sampled, and the error bars indicate the 
standard deviations. 



4.2 Reorientation and anisotropic compaction 

Rod-shaped clusters orienting to the rotational axis were 
observed in the simulations (see 3.3). Since the Brownian 
motion was not taken into account, the pioneer works by 



Jeffery on spheroids in shear flows may be referred to [48- 
50]. Spheroids can be considered as one of the simplest 
object representing elongated shapes. Jeffery analytically 
showed that, for the dilute limit, they should be in the pe- 
riodic orbits and have no tendency to set their axis in any 
particular direction under a simple shear flow. Besides, he 
expected that, since the dissipation of energy depends on 
their orbits, they would tend to adopt the orbital motion 
of the least dissipation of energy with additional elements 
in real suspensions. 

The trajectories of the principal axis for the smallest 
principal moment of inertia rii are plotted in Fig. 10 for 
two typical cases ending up as the rod-shaped (a) and 
round-shaped (b) clusters. Each of them shows the tra- 
jectory of one simulation with increasing shear rates. The 
color scale of the trajectories represents the changes of the 
shear rates. Though the orbits are not closed, one can flnd 
some similarity with the Jeffery orbits in the short time 
behavior (c.f. Figure 5.5 of ref. [50]). For the rod-shaped 
compaction (a), the orbit tends to converge in a narrow 
circuit around the north pole. It can be noticed that in 
the round-shaped compaction (b), the orbit is easier to be 
affected by the irregular structure of the cluster, in partic- 
ular when it crosses the x-y plane. Thus, if the orientation 
of clusters is changeable, the uniform compaction can be 
expected. 




Fig. 10. (Color online) The trajectories of the principal axis 
for the smallest principal moment of inertia ni are plotted on 
the unit sphere. The two examples are shown: (a) the rod- 
shaped compaction, which starts from {r,&/2-n) — (2.4,0.59) 
and ends to (2.7,0.18), and (b) the round-shaped compaction 
from (2.3,0.42) to (1.4,0.67). The north pole shows the rota- 
tional axis (y-axis). The color scale of the trajectories repre- 
sents the changes of the shear rates. 



Though the randomness of the cluster may lead to uni- 
form compaction, the formations of the rod-shaped clus- 
ters and their reorientations are not explained yet. An- 
other view point is the anisotropic compaction in shear 
flows. For clusters being restructured, the principal axis 
rii of the cluster is no longer fixed to the cluster, but de- 
pends on the structure at each instant. If the compaction is 
anisotropic, it may look as the reorientation of the princi- 
pal axis. In a shear flow, the drag forces acting on particles 
within rotating clusters increase with the distance from 
the rotational axis [ ' "]. So, the displacements of particles 
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depend on their positions as well. Since the compaction is 
caused by generations of cohesive bonds, one can expect 
that anisotropic compaction reduces distances of particles 
from the rotational axis. Thus, as long as the rotational 
axis is unchanged, the clusters tend to be compacted to 
elongated shapes. 



model and Dr. K. Ichiki for providing the simulator of Stoke- 
sian dynamics "Ryuon" and instructive advice, and Vincent 
Biirger for assistance with the manuscript. We acknowledge 
the financial support of the German Science Foundation (DFG 
priority program SPP 1273). 



4.3 Hydrodynamic effect 

In order to pronounce the hydrodynamic effect, the free- 
draining approximation (FDA) was compared with SD. 
First, the hydrodynamic effect is clearly seen in the hydro- 
dynamic stress for tenuous clusters, i.e. the critical shear 
rate with SD was much larger than the one with FDA: 



4.2 (see sect. 3.2). In low Reynolds num- 



p(SD) ,A(FDA) 
rc / rc 

ber flows, the disturbance of flows decays proportional to 
the inverse of the distance [49, 50], which results in the 
reduction of the drag forces acting on particles within iso- 
lated clusters [37]. Second, the difference between the two 
methods was also confirmed in the shape and orientation 
tendencies by the compaction (3.3.) The spatial distribu- 
tion of the drag force within clusters has the same sym- 
metry for FDA and SD [.37]. This is why some qualitative 
explanations for the formation of rod-shaped clusters are 
expected to be applicable for the simulation with not only 
SD but also FDA. However, as seen in Fig. 8, the result 
with FDA did not show clear tendency to form rod-shaped 
clusters. This result suggests that the hydrodynamic effect 
works as a kind of positive feedback for the anisotropic 
compaction in shear flows. 



5 Conclusion 

Restructuring of colloidal aggregates in shear flows has 
been investigated by coupling an interparticle contact 
model with Stokesian dynamics. We have introduced a 
method to reduce calculation cost for near-rigid behav- 
ior of aggregates in shear flows. That method has real- 
ized fluid-particle coupling simulations for a significant 
time period. The simulations with the stepwise increase of 
shear rate have demonstrated the reinforcement of clus- 
ters due to irreversible compaction with increase of the 
hydrodynamic stress. We expect that the observed con- 
solidation behavior induced by less-abrupt-application of 
flows is rather general. The introduced aspect ratio vs. ori- 
entation representation has also clarified the two types of 
compaction behaviors ending up in rod-shaped clusters 
orienting to the rotational axis and to round-shaped clus- 
ters. This anisotropic compaction can be a sort of hydro- 
dynamic effect, i.e. the enhancement of the tendency was 
seen in the comparison with the free-draining approxima- 
tion. Thus, the simulations with selected parameters for 
the contact model showed the characteristic behaviors of 
colloidal aggregates under flow conditions. 

The authors would like to thank Dr. Martine Meireles and Prof. 
Bernard Cabane for many valuable discussions on the contact 



A Descriptions by the equivalent ellipsoid 

In order to quantify the shape and orientation of random- 
structured clusters, we translate them to equivalent ellip- 
soids having the same principal moments of inertia [36]. 
The principal moments of inertia Ii, I2, I3 (/i < I2 ^ I3) 
are obtained by diagonalization of the moment of inertia 
tensor. By using them, the lengths of the semi-principal 
axes of the equivalent ellipsoid {a > b > c) are given as 
follows: 



5h_ 
2 



h - h 



N 



b = 



bh+h-h 
2 



N 



hh+h-h 
2 



N 



(31) 



The anisotropic shape of clusters is described by the as- 
pect ratio r of the equivalent ellipsoid: 



2a 



(32) 



In order to estimate the compaction of clusters, the effec- 
tive volume fraction t^eff is introduced as the ratio between 
the total volume of particles V-p = (47r/3)A^ and the vol- 
ume of the ellipsoid Vc = (47r/3) afoc, that is 



N 
abc 



(33) 
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